//============================================================================
// Name        : trjSelf.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "AtomsXYZ.h"
#include "HeaderTrj.h"
#include "FComms.h"
#include "Help.hpp"
#include "FstreamC.h"
#include "FstreamF.h"
#include "trjInput.h"


using namespace Topol_NS;
using namespace std;
unsigned int nnx=128,nny=128,nnz=128;
// profile variables

AtomsXYZ * atm;

int main(int argc, char ** argv)
{
	trj::trjInput MyInput(argc,argv);
	Parallel::comm=new Parallel::MyMPI(argc,argv);
	bool bDel=true;
	bool bWriteFab=false;
	bool bTestVol=false;
    bool bHyd=true;
	bool bIsrd=false;

	int narg=1;
	string filein,fileout="Micelles.out",filechg="Charges.dat",filepdb,filedomain,filedefdomain,fileirl,fileidb,fileclust="Cluster.out";
	string fileoutp1,fileoutp2,fileoutp3;

	string scratchdir;
	ifstream fpdb;
	ifstream ftest,fdefdomain;
	Fstream * finx=NULL;
	Fstream * fout_xtcx=NULL;
	ofstream * foutx_xyz=NULL;
	FstreamF * foutx=NULL;
	TopolPDB topPDB;

	ofstream fdomain;
	string comm0(argv[0]);
	size_t mypos=comm0.find_last_of("/")+1;
	size_t length=comm0.size()-mypos;
	string command=comm0.substr(mypos,length);
	string errmsg;
	string Usage="Usage:\t"+ command + "\n";
	vector<string> use=MyInput.getUsage();
	for(unsigned int n=0;n<use.size();n++)
		Usage+=use[n];
	Usage+="\n\t Default values in square brackets []\n";
	try{
		if(argc == 1) throw Usage;
		if(int m=MyInput.bTest().size()) {
			errmsg=" Command(s) not found: ";
			for(unsigned int n=0;n<m;n++)
				errmsg+=MyInput.bTest()[n]+"  ";
			errmsg+="\n"+Usage;
			throw errmsg;
			}
		}
	catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}

// Start input
	int nstart=0;
	int nend=-1;
	int nskip=1;
	// input read
	bool inputfile=false;
	try{
		if(!MyInput["-dcd"].empty()) {
			if(MyInput["-dcd"].size() < 2) throw string("\n Filename expected for " + MyInput["-dcd"][0] + " option \n");
			if(MyInput["-dcd"].size() > 2) throw string("\n More than one entry for " + MyInput["-dcd"][0] + " option \n");
			filein=MyInput["-dcd"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamF(filein,"rb");
			inputfile=true;
		}
		if(!MyInput["-xtc"].empty()) {
			if(MyInput["-xtc"].size() < 2) throw string("\n filename expected for " + MyInput["-xtc"][0] + " option\n ");
			if(MyInput["-xtc"].size() > 2) throw string("\n More than one entry for " + MyInput["-xtc"][0] + " option \n");
			filein=MyInput["-xtc"][1];
			ftest.open(filein.c_str(),ios::in);
			if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
			ftest.close();
			finx=new FstreamC(filein,"rb");
			inputfile=true;
		}
		if(!MyInput["-pdb"].empty()) {
			if(MyInput["-pdb"].size() < 2) throw string("\n filename expected for " + MyInput["-pdb"][0] + " option \n");
			if(MyInput["-pdb"].size() > 2) throw string("\n More than one entry for " + MyInput["-pdb"][0] + " option \n");
			filepdb=MyInput["-pdb"][1];
			fpdb.open(filepdb.c_str(),ios::in);
			if(!fpdb) throw string("\n Cannot open " + filepdb + "!!\n");
		}
		if(!MyInput["-o"].empty()) {
			if(MyInput["-o"].size() < 2) throw string("\n filename expected for " + MyInput["-o"][0] + " option \n");
			if(MyInput["-o"].size() > 2) throw string("\n More than one entry for " + MyInput["-o"][0] + " option \n");
			fileout=MyInput["-o"][1];
		}
		if(!MyInput["-b"].empty()) {
			if(MyInput["-b"].size() != 2) throw string(" Number of first frame needed for " + MyInput["-b"][0] + " option ");
			stringstream(MyInput["-b"][1])>> nstart;
		}
		if(!MyInput["-e"].empty()) {
			if(MyInput["-e"].size() != 2) throw string(" Number of last frame needed for " + MyInput["-e"][0] + " option ");
			stringstream(MyInput["-e"][1])>> nend;
		}
		if(!MyInput["-skip"].empty()) {
			if(MyInput["-skip"].size() != 2) throw string(" Number of skipped frames needed for " + MyInput["-skip"][0] + " option ");
			stringstream(MyInput["-skip"][1])>> nskip;
		}
		if(!MyInput["-help"].empty()) {
			if(MyInput["-help"].size() != 1) throw string(" No argument for " + MyInput["-help"][0] + " option ");
				throw Help;
		}
		if(!fpdb.is_open()) throw string("\n Need to provide .pdb file !!\n");

	}catch(const string & s){
		cout << s << endl;
		Parallel::comm->Barrier();
		Parallel::comm->~MyMPI();
		return 0;
	}
	if(fileout.substr(fileout.find_first_of(".")+1).find("txt") != string::npos)
		foutx_xyz=new ofstream(fileout.c_str(),ios::out);
	else
		foutx=new FstreamF(fileout.c_str(),"wb");

	atm=new AtomsXYZ;


	Topol::fDomainPDB=new std::ofstream ;
	Topol::fDomainPDB=&fdomain;


	Fstream & fin=*finx;
	Fstream & fout_xtc=*fout_xtcx;
	FstreamF & fout=*foutx;
	ofstream & fout_xyz=*foutx_xyz;


	try{
		if(!inputfile) throw string("\nNo input file given will operate on the .pdb file\n");
	}
	catch(const string & s ){
		cout << s << endl;
		Parallel::comm->Barrier();
	}

	// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}
	// Define radii to use

	TopolMic MyTop;


	// def domain
	if(Topol::bDefDomain) {
		string Dom0;
		int begin0,end0;
		vector<int> tmp(2);
		while(fdefdomain >> Dom0 >> tmp[0] >> tmp[1]) {
			MyTop.AddDomain(Dom0,tmp);
		}
	}

	topPDB(data);
	MyTop(topPDB,false);
//	atm->setrd(MyTop);

	typedef map<string,vector<vector<int> > > mappa;
	mappa & Def=MyTop.getDef();
	vector<string> CSelect;
// Micelle stuff

	int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();


	atm->setDim(natoms);
	ios::streampos len;
	HeaderTrj header;

// Read header of dcd file
	if(finx) {
		fin.seekg(0,"end");
		len=fin.tellg();
		fin.seekg(0,"beg");
		fin>>header;
		try{if(!header.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not not match!");}
		catch(const string & s){cout << s<<endl;exit(-1);}
	}
	atm->initPDB(topPDB);


// Read trajectory and compute Micelle data

    Parallel::comm->Barrier();
    AtomsXYZ b=*(static_cast<AtomsXYZ *> (atm));
	int ntime=0;
	int nnStep=0;
	while(fin.tellg() < len){
		if(ntime < nstart) {fin+=*atm;ntime++;continue;}
		if(nend != -1 && ntime > nend) break;
		fin >> *atm;
		if(!(ntime%nskip)) {
			cout << fixed << setw(5) << "----> Step No. " << setw(5) << right << fixed << ntime << " -- Time = " << atm->getTime() <<" ps \n";
			if(foutx_xyz){
				fout_xyz << *atm;
			} else if(foutx){
				fout<< *atm;
			}
		}
		Parallel::comm->Barrier();
		ntime++;
		nnStep++;
	};

	Parallel::comm->Barrier();
	Parallel::comm->~MyMPI();
	cout << "\nProgram completed: Output data written to " + fileout << "\n\n";

	if(foutx_xyz){
		fout_xyz.close();
	} else if(foutx){
		fout.getfout().close();
	}
    return 0;
}   //end of main function
